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PROGRESS  IN  THE  MODELING  OF  THE  SHOCK  RESPONSE 
AND  MITIGATON  OF  THICK  COMPOSITE  SHELLS 


1.  INTRODUCTION 

Future  naval  vessals  are  expected  to  utilize  increasing  iimounts  of  thick  polymer  matrix  com¬ 
posite  materials  in  their  structure  due  to  significant  system  performance  enhancements  achievable 
with  tliese  stiff,  strong,  lightweight  materials.  In  order  to  apply  these  materials  to  submarine  and 
surface  ship  structures,  their  behavior  under  highly  transient  shock  loadings  must  be  better  under¬ 
stood.  Current  analysis  methods  for  impact  and  underwater  shock  response  were  developed  pri¬ 
marily  to  predict  the  behavior  of  ductile,  homogeneous  metal  structures,  which  under  severe  loads 
deform  plastically  due  to  slip  along  shear  planes.  Composites  are  very  heterogeneous  materials, 
combining  high  performance  fibers  in  a  viscoelastic  matrix.  Also,  composites  tend  to  be  highly 
dispersive  to  propagating  waves  and  deform  nonlineardy  under  severe  load  through  the  develop¬ 
ment  of  networks  of  micro-cracks. 

This  report  is  a  continuation  of  our  efforts  [1-5]  to  develop  a  methodology  for  predicting  the 
response  of  thick  composite  materials  subjected  to  multi -dimensional  shock  loadings.  In  develop¬ 
ing  the  numerical  aspects  of  this  work,  the  explicit  transient  finite  element  programs  PR0NT02D 
[8]  and  ABAQUS/EXPLICIT[9]  are  being  employed  to  construct  two  dimensional  (2D)  and  three 
dimensional  (3D)  capabilities.  In  addition,  the  explicit  finite  difference  program  WONDY[7]  is 
being  used  to  study  one  dimensional  (ID)  behavior  and  to  develop  advanced  material  models  re¬ 
lating  to  evolving  damage,  dispersion  and  viscoelasticity.  Appendix  A  contains  a  brief  description 
of  the  reduced  ID  damage  model  with  nonlinear  compression,  which  has  been  programmed  into 
WONDY.  This  constitutive  model  is  based  upon  the  2D  continuum  damage  model  developed  in 
[1-3). 

In  section  2,  the  inclusion  of  dispersion  effects  in  the  ID  damage  in  the  damage  model  is  dis¬ 
cussed.  The  use  of  higher  derivatives  in  the  constitutive  modeling  is  explored  as  well  as  predictions 
for  composite  plates  subjected  to  underwater  shock.  The  effects  of  viscoelasticity  in  the  matrix  of 
the  composite  plate  and  its  effect  on  dispersion  is  also  addressed. 

The  application  of  the  2D  continuum  damage  theory  is  discussed  in  section  3.  A  comparison 
is  made  'a  PR0NT02D  models  to  ID  damage  results  from  WONDY  predictions  for  tlie  impact 
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of  a  plexiglass  flyer  and  graphite/peek  plate  -  see  [6]  for  the  expcriniei\tal  results.  This  is  followed 
by  a  coarse  model  of  2D  impact  using  the  PR0NT02D  program. 

In  section  4,  a  3D  continuum  damage  theory  is  developed  for  thick  laminated  composites.  The 
3D  formulation  in  part  is  an  extension  of  the  2D  transversely  isotropic  damage  theory  [1-3]. 
However,  three  dimensional  considerations  as  well  as  the  inefficiencies  of  trying  to  model  indi¬ 
vidual  plies  in  a  thick  composite  require  a  slightly  different  approach.  In  this  section,  a  3D  formu¬ 
lation  is  developed  which  applies  the  damage  directly  to  the  stresses  as  in  [10-12],  rather  than  the 
compliances  as  is  the  case  in  the  2D  theory. 

ABAQUS/EXPLICIT  [9]  will  in  the  future  be  the  major  developmental  code  for  2D  and  3D 
finite  element  analysis  replacing  PRONTO.  Appendix  B  is  a  listing  of  the  ABAQUS/EXPLICIT 
user  written  subroutine  for  2D  continuum  damage,  and  is  similar  to  the  2D  damage  model  pro¬ 
grammed  into  PR0NT02D. 

Finally  in  sections  5  and  6,  a  summary  is  given  followed  by  some  thoughts  on  future  direc¬ 
tions. 
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2.  DISPERSION  EFFECTS  IN  A  ID  DAMAGE  MODEL 


2.1  Background  and  general  remarks 


Composite  materials,  especially  laminated  plates  and  shells,  comprise  a  very  heterogeneous 
media.  The  dynamic  response  of  such  materials  can  be  broadly  classified  into  one  of  two  group¬ 
ings  [13].  If  the  wavelength  of  the  loading  and  the  response  of  the  material  is  very  long  compared 
with  the  scale  of  the  inhomogeneity,  then  the  material  response  is  governed  by  effective  properties 
of  the  equivalent  homogenized  media.  However,  if  the  composite  structure  is  subjected  to  shock 
loadings  such  as  UNDEX  or  impact,  then  the  wavelengths  of  the  loading  and  the  response  of  the 
structure  are  much  shorter.  If  this  is  the  case,  the  characteristic  dimensions  of  the  hetergeneous  me¬ 
dia  become  much  more  important  and  the  dynamics  is  greatly  complicated.  The  interfaces  be¬ 
tween  the  material  phases  cause  wave  reflection  and  refraction.  The  energy  is  thus  spread  or 
“dispersed”  over  many  wavelengths.  Mathematically  this  implies  that  the  phase  speed  c  and  the 
frequency  o)  are  a  function  of  the  wavelength  k  [13-14]  or: 


C((D) 


m(k) 

k 


(2Aa) 


In  a  no.i-dispersive  media,  the  phase  speed  equals  the  wave  speed  and  is  a  constant.  The  frequency 
is  then  expressed  as: 

(d  =  CQk  (,12b) 

Dispersion  in  composite  materials  can  have  two  sources.  We  have  already  discussed  the  role 
that  geometry  changes  can  play  in  a  heterogeneous  media.  Another  form  of  this  phenomcnia  has 
as  its  origins  the  viscoelastic  nature  of  the  matrix  of  the  composite.  Separating  the  geometric  and 
the  viscoelastic  effects  is,  however,  a  very  difficult  undertaking,  in  [15]  Sutherland  discusses  this 
very  problem. 

For  our  purposes,  however,  it  is  not  necessary  in  the  numerical  methods  to  completely  separate 
these  two  sources  of  dispersion.  Rather  what  is  needed  is  an  approach  that  can  accurately  take  into 
account  dispersion  in  general  without  requiring  an  excessive  number  of  degrees  of  freedom  and 

associated  computer  costs. 
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Before  describing  the  specific  details  of  our  approach,  we  note  that  as  discussed  in  [14],  the 
phenomenon  of  dispersion  can  be  introduced  directly  into  the  numerical  schemes  such  as  finite  dif¬ 
ferences  and  finite  element  methods.  For  homogeneous  media,  numerical  dispersion  is  not  desired 
since  it  represents  additional  error  and  inaccuracy  in  the  response  of  the  system  being  modeled. 
However  as  indicated  in  [14,16-17],  unwanted  numerical  dispersion  is  usually  present  in  computer 
models. 


2.2  Development  of  an  approach  and  application 


For  a  heterogeneous  media,  the  problem  encountered  involves  tlie  inclusion  of  the  dispersive 
effects  without  the  need  to  explicitly  model  all  the  fine  details  of  the  materials,  which  would  re¬ 
quire  a  very  large  number  of  degrees  of  freedom.  One  approach  is  to  introduce  higher  derivatives 
into  the  governing  equations  [17-18].  In  [17],  longitudinal  waves  in  a  plate  are  studied.  The  clas¬ 
sical  wave  equation  is: 


_  20 


(2.2) 


where  u  is  the  longitudinal  displacement  and  Cq  is  the  wave  speed  (non-dis|>ersive).  To  include 
dispersion  effects,  a  fouilh  order  partial  differential  equation  is  substituted  in  place  of  (2.2)  that  has 
the  form; 

+  =0  (2.3) 

dr  dx^ 

where  ot,  P,  y  and  k  arc  constants  chosen  to  satisfy  the  governing  dispersion  equation. 


An  approach  similar  to  [17]  lias  been  incorporated  into  the  one  dimensional  program  WONDY 
in  an  attempt  to  include  dispersion  effects  due  to  the  heterogeneity  of  the  composite  material 
[4,19].  For  a  ID  case,  the  stress  strain  relationship  can  be  written  as: 

C  =  C,,€  (2.4) 

Next  higher  derivatives  with  respect  to  time  arc  introduced  into  eq.  (2.4)  in  the  form: 


Y,C|iC)  +  f^  +  A](0|-T2C,,F.)  =0 

vN  '^2/ 


(2.5) 
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where  A,j,  Xj,  Yj  and  are  empirical  constants  to  be  determined  by  matching  with  experimental 
results.  Equation  (2.5)  is  a  relatively  straight  forward  constitutive  relation  tind  its  numerical  inte¬ 
gration  within  an  explicit  program  such  as  WONDY  would  appear  to  be  routine.  In  fact,  the  inte 
gration  using  Runge-Kutta  methods  of  eq.  (2.5)  proved  to  be  quite  challenging,  requiring  even  finer 
time  steps  than  needed  by  the  central  difference  operator  used  in  WONDY  to  integrate  the  equa¬ 
tions  of  equilibrium.  The  use  of  explicit  numerical  integration  for  the  constitutive  relation  eq.  (2.5) 
thus  proved  to  be  too  expensive. 

Another  approach  to  the  integration  of  eq.  (2.5)  is  to  express  it  in  a  convolution  product  form 
such  as: 


cr=(G®e)  (2.6) 

where  <8>  represents  the  convolution.  Toward  this  end,  one  can  take  the  Laplace  transform  of 
cq.  (2.5)  assuming  zero  initial  conditions  [19].  Doing  this  and  then  performing  an  inverse  Laplace 
transformation  produces; 


a(t\  = 


(">  1\ 


where  H{t)  is  the  Heaviside  function,  and  G(t-t')  is  a  function  of  I,  t\  y,,  Y2,  X,  and  X,  •  see 
[19],  Equation  (2.7)  has  the  general  form  of  the  hereditary  integrals  produced  in  viscoelasticity 
theoi'y  [20].  The  stress  at  time  j  can  be  expressed  by; 


^  + 1 )  =  72^  ( 1 )  +  (x;  ( 27,  -  Y2  -  1 )  sin 

«+  I  -L 
>  =  i 


( 1  -  Y  )  cos 


V  ^2  // 


j^(2Y, -Yj- l)cos|^^j  (1-Y2)sin 


\  '^2  J 


w-f  1  ± 

X 

y  =  1 


(2.8) 


where 
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(2.9fl) 


^)dt  {2.9b) 

The  variables  and  /,  in  eqs.  (2.9)  can  be  explicitly  integrated.  The  stress  at  time  -n .  O  (/„  + , ) 
then  becomes  a  matter  of  summing  the  hereditary  integrals  over  the  previous  time  steps.  This  is 
much  more  computationally  efficient  than  the  numerical  integration  of  the  differential  equation 
(2.5). 

A  great  many  WONDY  computer  nins  simulating  the  plexiglass  flyer  and  graphite/peek  plate 
impact  experiments  in  [6]  were  performed.  Equations  (2.5),  (2.7)  -  (2.9)  were  used  to  model  the 
dispersion  effects  in  the  graphite/peek  plates.  Several  of  the  experiments  were  modeled  in  an  at¬ 
tempt  to  determine  the  enipircal  constants  Xj,  Tj  Y2  sufficient  accuracy  for  graphite/ 
peek  material.  An  important  criteria  was  the  matching  as  close  as  possible  to  the  experiments  of 
the  plate  backface  particle  velocity.  Curve  fitting  to  separate  ex|)erimeuts  produced  too  much  scat¬ 
ter  in  the  enipimnl  constanis.  This  indicaieti  that  the  origin,.!  disjxjrsion  assumption  in  the  form  of 
eq.  (2.5)  was  not  fully  satisfactory. 

In  an  attempt  to  improve  the  material  model,  the  following  viscoelastic  type  of  constitutive  law 
was  postulated  to  account  for  dispersion: 

t 

a(0  =  //(r)jG(r--r')e(/)r/r  (2.10) 

0 

where 

-0.05  (-)  , 

O  (t)  =  yCjj  -1-  ( 1  “  7)  Cjj?  cos  ( ^)  (2. !  1) 

This  constiutive  model  yielded  slightly  better  results  than  the  previou.s  one  -  eqs.  (2.5)  and  (2.7)  - 
(2.9).  In  Table  2.1  the  modeling  of  impact  experiment  no.  204  [6],  which  consisted  of  a  plexiglass 
flyer  and  a  graphite/peek  plate,  with  WONDY  is  described.  See  Figure  3.4  for  a  description  of  the 
impact  experiments.  Figure  2.1  indicates  the  results.  The  predicted  and  experimental  backface 

A 


0+1  -hlLJl 

=  Ez  ^  j  e  cos( 

P*' =  Ei  ^  Jc  '  sin( 


particle  velocity  of  the  plate  is  shown  in  Figure  2.1a.  In  general,  the  agreement  is  very  good,  with 
some  deterioration  later  in  the  simulation  ns  might  be  expected.  Figure  2.1b  describes  the  predict¬ 
ed  Vj  (thru  thickness)  damage  in  tiie  plate.  A  value  of  1  indicates  total  delamination. 

Overall,  the  viscoelastic/dispcrsion  constitutive  law  represented  by  cqs.  (2.10)  and  (2.1 1)  of¬ 
fered  some  improvement  over  the  original  disi)ersion  model  expressed  by  eq.  (2.5).  Cur\'e  fitting 
to  separate  ext>eiiments  produced  less  scatter,  but  in  general  the  results  weie  not  fully  saiisfactoi  y 
-  the  constitutive  model  was  just  not  robust  enough. 

In  reference  [4],  the  authors  programmed  the  viscoelastic  constiutive  law  discussed  in  [15] 
into  WONDY,  and  modeled  some  of  the  impact  exi)criments  in  [6]  which  consisted  of  plexiglass 
flyers  and  graphite/peek  plates.  Their  results  appear  to  be  very  good  with  less  scatter  among  the 
numerical  results  for  different  flyer/plate  thicknesses,  initial  velocities,  etc.  than  was  obtained  in 
this  report.  This  approach  [4]  will  probably  be  incorporated  in  future  2D  and  3D  continuum  dam¬ 
age  work. 
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♦  3.45  mm  PMMA  flyer  at  93  in/s  onto  25.83  mm  Gr/PEEK 
target  (both  flyer  and  target  102  mm  diameter). 

•  Mesh  Size— 215  zones  in  flyer  and  1500  zones  in  target. 

♦  Geometric  dispersion-dissipation  constitutive  law 

o  =  H(t) |oG(t  -  t‘)e(x,  t')dt  with  G{t)  =  yC,,  +  (1  -  y)  C,,o‘”*‘''’'cos  (t/t). 

*  Coupling  between  dispersion  and  compressive  nonlinearity 

Y-  {Yo  +  6E)/(1+8E)  where  8E  =  Ax 

•  Dispersion  properties  x=0.01  psec,  yssO.? 

♦  Coupling  constant  A= 150.0 

*  Damage  properties  rii=2,5  psec,  ni=s2,  a^oss  /O  MPa  where 

-V?)* 


Table  2.1  Modeling  of  impact  experiment  no.  204  [6],  plexiglass  flyer  and 
graphite/pcck  plate. 
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Particle  Velocity  m/sec 


FlGUItE  2.1 


SPECIMEN  NO.  204-PREDICTIONS  AND  COMPARISON  WITH  EXPERIMENT 
(PMMA  AT  93  M/S  ONTO  GRAPHITE/PEEK) 


n.  Backfncc  Purticic  Vclocity-Prcdiclion  b.  Pk'cdictcd  Vj  Daiiingc 

Compared  wlUi  ExpcrimciUal  Data 


lU  U  M  U  lt>  0.2  0.4  0.6  0.0  1 

Time  After  Impact  pscc  Dimensionless  Plate  Thickness 
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23  Predictions  for  composite  plates  subjected  to  undervt^ater  shock 

In  this  section,  the  WONDY  program  is  employed  to  predict  the  response  of  composite  graph- 
ite/peek  plates  from  2  to  8  inches  which  are  subjected  to  simulated  near-field  underwater  explo¬ 
sions.  Dispersion  and  damage  are  included  in  the  material  models  used  in  the  WONDY  computer 
runs.  The  intent  of  these  analyses  is  to  assess  the  extent  of  spallation  damage  for  various  charge 
weights  and  ranges,  and  to  determine  the  approximate  boundaries  of  incipient  and  complete  spal¬ 
lation.  The  history  of  a  TNT  explosive  impulse  passing  through  a  fixed  spatial  point  can  be  ex¬ 
pressed  as  [21]; 


P  =  ®  (2.10) 

where  pQ  is  the  pressure  asid  0  the  exponential  decay  constant.  According  to  Aaron’s  law  and 
0  can  be  determined  using: 

y^l/3 

Po  =  0.519(^)  (2.11) 

y,l/3  -022 

0  -  0.925  (-51^) 

A 

where  Pq  is  given  in  kilobars,  0  is  in  seconds,  R  is  the  range  in  meters  (m),  and  W  is  the  charge 
mass  in  kilograms  (kg).  The  values  of  H'  in  eq.  (2.1 1)  will  range  from  5  kg  to  30  kg,  R  will  vary 
from  .5  to  l.S  m  (very  close  in),  and  the  plates  as  previously  stated  will  vary  from  2  to  8  inches  or 
.0508  to  .1016  meters. 

For  seawater,  the  following  material  properties  (SI  units)are  employed  in  the  modeling: 

p  =  1000^ 

=  1500-  (2.12) 

K  -  1-75 

where  p  is  the  density  and  Cq  is  the  bulk  sound  speed.  The  variable  s  comes  from  a  linear  shock 
velocity  to  particle  velocity  relationship  of  the  form  [7,8]: 
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f/,  =  (2.13) 

The  material  constants  used  for  composite  plates,  which  are  assumed  to  be  composed  of  giaph- 
ite  fibers  and  a  peek  matrix  (see  [6]  and  Appendix  A)  are  as  follows: 

bulk  properties 

p  =  1579^  (2.14a) 

m 

Co  =  3000- 
^  s 

£,j  =  69.0  X  10^ Pa  (thru  thickness) 

£22  =  13-45  X  lO^Pa 
v,2  =  0.(K,  Vj3  =  0.3 

damage  parameters  (see  [1-3]) 

n.  =  1.0x10“*^,  =  1.0,  ~  lOrAO^Pa  (2.14b) 

and  the  dispersion  parameters  (using  the  early  dispersion  model  -  see  eq.  (2.5)) 

X,  =  -  ,  Xj  =  5.0X  lO"*^.  Yi  =  -  .and  =  0.81  (2.14c) 

Figure  2.2  describes  the  WONDY  model  employed  for  the  2  inch  thick  plate  and  the  fluid  me¬ 
dia.  We  note  that  the  model  for  the  fluid  domain  is  only  one  half  the  thickness  of  the  plate.  In  gen¬ 
eral  it  is  desirable  to  include  as  little  of  the  fluid  as  possible  since  the  pressure  wave,  whicii  is 
generated  from  a  standoff  position  ranging  from  .5  to  1.5  meters  in  these  analyses,  must  pass 
through  the  fluid  on  its  way  to  the  solid.  The  bulk  sound  speed  of  the  fluid  is  one  half  that  of 
the  composite  Cq.  Thus  the  time  of  flight  (=// c)  is  roughly  the  same  for  both  the  fluid  and  the 
composite  media.  In  addition  in  these  analyses,  it  was  not  necessary  to  employ  silent  boundary  con¬ 
ditions  [7,8]  at  point  A  of  the  fluid  boundary.  The  analyses  were  terminated  (the  damage  was  done 
by  the  reflected  tensile  wave  from  the  backface,  point  B,  of  the  composite  plate)  before  any  spuri¬ 
ous  reflections  from  point  A  could  reach  the  plate. 
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One  other  important  point  to  be  brought  up  with  respect  to  Figure  2.2  concerns  the  number  of 
finite  difference  zones  used  in  the  fluid  and  the  composite  plate.  For  the  solid,  1000  zones  were 
typically  used  to  capture  the  damage  and  to  accurately  model  the  shock  wave  as  it  moved  through 
the  plate.  This  dicated  the  using  of  500  zones  in  the  fluid  media.  Employing  zoning  too  coarse 
(with  respect  to  the  solid)  in  the  fluid  would  result  in  the  shock  wave  not  being  accurately  propa¬ 
gated  through  the  fluid.  The  selection  of  the  time  step  5r  is  based  upon  the  smallest  l/c  (/  is  the 
length  of  a  zone)  in  the  model,  which  for  the  model  described  in  Figure  2.2  will  be  determined  by 
the  composite  plate.  Zoning  too  coarse  in  the  fluid  would  not  allow  the  shock  wave  to  get  across 
a  zone  in  one  time  step  and  cause  a  diffusion  of  the  wave. 

In  general,  the  use  of  1500  zones  in  the  finite  difference  meshes  for  this  fluid-structure  inter¬ 
action  problem  would  seem  to  be  excessive.  Based  on  the  discussion  in  section  3.2  this  may  very 
well  be  the  case.  However,  this  point  is  not  clear  at  this  time  and  more  modeling  needs  to  be  per¬ 
formed  to  determine  minimum  levels  of  discretization  required  to  accurately  model  both  the  shock 
wave  and  the  overall  damage  profile. 

Finally  in  Figure  2.3,  the  results  are  shown  for  2, 4  and  8  inch  thich  graphite/peek  plates  sub¬ 
jected  to  underwater  explosions.  Based  on  these  predictions,  it  appears  that  spallation  is  not  a  real 
threat  except  for  very  close  in  explosions. 
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P:=Poe 


1  inch  (.0254  m) 


2  inch  (.O508m) 


FLUID  MEDIA  (500  zones) 


fluid-structure  interface 


SOUD  MEDIA 


GRAPHITE/PEEK  PLATE 


(1  (XX)  zones) 


B  (baciface  of  plate) 


FIGURE  2.2  WONDY  model  for  composite  plate  subjected  to  underwater  shock. 
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FIGURE  2  J 


EXPLOSION-INDUCED  SPALLATION  DAMAGE 
CALCULATIONS  IN  A  CHARGE  WEIGHT  VS  CHARGE  RANGE  SPACE 
FOR  TNT  IN  WATER  OVER  GRAPHITE/PEEK 
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3.  FURTHER  APPLICATIONS  OF  THE  2D  TRANSVERSELY  ISOTROPIC 
CONTINUUM  DAMAGE  MODEL 

3.1  Background 

Details  of  the  2D  transversely  isotropic  continuum  damage  model  developed  at  NRL  arc  con¬ 
tained  in  references  [1-3].  Also,  section  4.1  of  this  report  briefly  discusses  some  of  the  salient  fea¬ 
tures  of  the  2D  theory  in  developing  a  viable  3D  continuum  damage  approach. 

Basically  our  efforts  to  run  detailed  2D  continuum  damage  models  with  PR0NT02D  have 
been  hindered  recently.  We  have  begun  the  transition  of  the  2D  and  future  3D  damage  woric  away 
from  PRONTO  to  the  program  ABAQUS/EXPLICrr  [9].  This  latter  explicit  finite  element  pro¬ 
gram  has  been  developed  by  the  architects  of  PRONTO.  ABAQUS/EXPLICIT  thus  represents 
second  generation  software  that  possesses  good  pre-  and  post-processing  capabilities,  and  is  well 
documented  and  supported. 

Another  event  that  has  hindered  recent  efforts  to  run  detailed  2D  damage  models  has  been  the 
retirement  of  the  NRL  CRAY-XMP,  which  was  a  COS  machine.  PR0NT02D  was  run  in  the  past 
on  this  COS  machine  for  larger  jobs.  A  CRAY-EL  is  now  available  at  NRL  with  the  UNICOS  op¬ 
erating  system,  but  our  development  version  of  PR0NT02D  could  not  be  converted  due  to  the 
COS  libraries  which  are  linked  to  the  source  code.  Rather  than  trying  to  update  a  new  UNICOS 
version  of  PRONT02D,  a  decision  was  made  to  go  completely  to  ABAQUS/EXPLICIT  since  it 
would  be  employed  to  develop  3D  continuum  damage  capabilities.  Appendix  B  is  a  listing  of  a 
user  written  subroutine  for  the  current  2D  damage  theory  being  implemented  into  ABAQUS/EX¬ 
PLICIT. 

In  the  remainder  of  this  section,  we  discuss  PR0NT02D  analyses  which  were  run  on  a  Micro- 
Vax,  and  are  therefore  limited  in  the  fineness  of  the  meshes  and  the  number  of  elements  allowed. 
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3.2  PR0NT02D  comparisons  to  experimental  and  WONDY  results 

The  intent  of  this  section  is  to  apply  the  2D  continuum  damage  theory  of  [1-3],  which  has  been 
programmed  into  PR0NT02D,  to  a  one  dimensional  (ID)  impact  problem  that  has  been  studied 
experimentally  in  [6]  and  with  WONDY  *  see  Figure  2.1  for  results.  As  discussed  in  section  2.1 
and  Appendix  A,  the  WONDY  program  is  applicable  to  one  dimensional  situations,  and  contains 
both  a  reduced  ID  damage  model  and  a  capability  to  account  for  dispersion/viscoelastic  effects. 

We  address  two  questions  in  this  section.  The  first  concerns  the  effects  of  dispersion/vis¬ 
coelasticity  in  the  analyses,  which  is  not  explicitly  included  in  the  PRONT02D  model  but  is  pro¬ 
grammed  into  WONDY.  The  second  question  relates  to  the  number  of  finite  elements  neededwith 
PR0NT02D  (or  ABAQUS/EXPLICIT)  in  the  thru  thickness  direction  to  accurately  predict  the 
damage  profile  and  the  shock  wave.  Previous  WONDY  finite  diffemce  runs  had  seemed  to  indi¬ 
cate  that  hundreds  of  finite  difference  zones  were  necessary  to  accurately  simulate  the  shock  wave 
thru  the  thickness  of  the  plate.  For  instance,  the  WONDY  results  of  Figure  2.1  employed  1500 
zones  thru  the  composite  plate.  For  2D  and  especially  3D  problems,  this  would  tranlate  to  hun¬ 
dreds  of  finite  elements  thru  the  thickness,  and  is  completely  impractical.  Due  to  aspect  ratio  con¬ 
siderations,  tens  of  thousands  of  elements  would  be  required  to  model  just  a  2D  plate  an'*  flyer. 

With  this  in  mind,  we  therefore  focus  upon  the  PRONT02D  modeling  of  impact  experiment 
no.  204  in  reference  [6].  Table  2.1  and  Figure  2.1  describe  the  WONDY  modeling  and  results 
along  with  the  experimental  results  for  this  impact  problem,  which  consists  of  a  plexiglass(PMM  A) 
flyer  and  a  graphite/peek  plate  of  the  same  diameter.  Because  both  the  flyer  and  the  plate  aie  the 
same  diameter,  the  early  time  response  of  the  center  of  the  plate  can  be  approximated  using  a  one 
dimensional  model.  Once  bending  and  tranverse  shear  waves  reach  the  center  portion  of  the  plate, 
this  assumption  is  of  course  no  longer  valid.  Table  3.1  indicates  the  material  constants  for  the 
PMMA  flyer  along  with  the  material  and  damage  variables  employed  for  the  graphite/peek  plate. 
(See  also  Figure  3.5,  w'hich  indicates  the  material  constants  used  in  a  typical  PR0N1'02D  run.  For 
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a  further  explaination  of  the  variables  used  in  the  2D  continuum  damage  model  programmed  into 
PRONT02D,  refer  to  [1-3,8]).  Two  analyses  are  performed.  Model  A  employs  tOO  elements  thru 
the  plexiglass  flyer  and  400  elements  thru  the  graphite/peek  plate.  Model  B  is  coarser  and  uses  a 
quarter  the  number  of  elements,  or  25  thru  the  flyer  and  100  tiuu  the  composite  plate.  Also  note 
that  model  B  represents  a  viable  nK^sh  spacing  in  the  thickness  direction  that  could  be  extended  for 
a  troc  2D  model  of  this  problem,  including  bending  and  shear  effects  that  would  arise  later  in  the 
analysis.  (See  section  3.3  for  a  coarse  2D  simulation  of  this  problem  with  PR0NT02D.) 

The  PR0NT02D  damage  predictions  for  models  A  and  B  are  indicated  in  Figures  3.2.1  and 
3.2.2,  respectively.  Delamination  damage  in  the  plate  is  indicated  in  both  figures.  Comparing  Fig¬ 
ure  2.1  to  these  results,  we  see  that  model  B  (coarse  mesh)  produced  damage  predictions  reason¬ 
ably  close  to  the  veiy  refined  WONDY  model  (1500  zones).  Model  A  with  400  elements  thra  the 
thickness  of  the  plate  was  more  ragged  in  the  description  of  the  total  damage  (V]=:l)  than  model  B. 

In  general,  the  delamination  results  of  models  A  and  B  indicate  that  the  coarser  mesh  B  com- 
pai^d  better  to  the  very  refined  WONDY  model.  At  this  point,  it  is  thought  that  the  coarser  finite 
element  (model  B)  may  have  compensated  somewhat  for  the  dispersion/viscoelasticity  effects 
which  are  not  included  in  PR0NT02D.  (See  also  the  comments  concerning  Figures  3.3.2  and 
3.3.3).  This  is  only  speculation  using  at  this  point  in  time  the  WONDY  results  as  the  standard.  In 
the  future,  comparison  will  be  made  to  the  data  from  ultrasound  tests  [6],  and  the  forthcoming  dis¬ 
section  of  the  plates. 

In  Figure  3.3.1,  the  predicted  backface  velocity  of  the  plate  from  the  coarse  mesh,  model  B,  is 
indicated.  Comparing  to  the  experimental  results  which  are  shown  in  Figure  2.1,  one  sees  a  good 
correlation.  The  peaks  denoted  by  the  points  A,  B,  C,  and  D  in  Figure  3.3.1  match  up  very  well 
with  the  experimental  results.  The  width  of  the  initial  shock  front,  however,  is  narrower  for  the 
coarse  PR0NT02D  model  and  model  B  does  not  “shock  up”  as  quickly  as  the  experimental  curve. 
Points  E,  F  and  G  in  Figure  3.3.1  also  do  not  match  up  well  with  Figure  2.1.  The  predicted  back- 
face  velocity  is  somewhat  out  of  phase.  Overall,  however,  the  results  of  Figure  3.3.1  nre  very  en¬ 
couraging  for  our  2D  (and  3D  -  section  4)  continuum  damage  theories.  The  inclusion  of  viscoelas- 
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tic/dispersion  capabilities  should  greatly  assist  in  the  closing  of  the  gap  between  the  experimental 
and  predicted  results,  and  offer  a  real  hope  in  the  modeling  of  2D  and  3D  impact  problems. 

In  Figure  3.3.2,  the  predicted  backface  velocity  of  the  composite  plate  is  shown  for  the  fine 
mesh,  model  a  (400 elements  thru  the  thickness).  Figure  3.3.3  displays  the  results  from  both  model 
A  and  model  B  (coarse  mesh).  Comparing  to  the  experimental  results  in  Figure  2.1,  we  see  that 
the  fine  mesh  A  **shocks  up”  quicker  than  mesh  B.  The  initial  shock  front  is  also  wider,  and  thus 
closer  to  the  experimental  results.  However,  pa.st  approximately  1 1  microseconds,  the  fine  mesh 
predictions  deteriorate  and  eventually  around  18  microseconds  severe  oscillations  are  evident.  (It 
is  noted  that  with  400  elements  thru  the  thickness  of  the  plate,  mesh  A  requires  time  steps  approx¬ 
imately  1/4  those  of  mesh  B  to  maintain  stability  based  on  the  Courant  criteria  [8,9,23].  This  can 
lead  to  accumulated  error  later  in  the  time  history,  but  should  not  cause  such  instability  in  general.) 

Artificial  viscosity  [8]  is  contained  in  PR0NT02D  to  prevent  high  velocity  gradients  from 
collapsing  an  element  before  it  can  respond,  and  to  help  prevent  oscillations  or  “ringing”  such  as 
seen  in  Figures  3.3.2  and  3.3.3.  The  following  values  for  aiticfical  viscosity  were  used  in  obtaining 
the  results  shown  in  Figure  3.3.2: 

h,  =0.1  and  b2  =  2. 

where  b  j  is  the  linear  bulk  viscosity  coefficient  and  is  the  quadratic  viscosity  coefficient.  As 
discussed  in  [8],  helps  to  dissipate  truncation  frequency  oscillations  while  b2  is  designed  to 
smear  a  shock  wave  across  several  elements.  In  an  attempt  to  decrease  the  oscillations,  h^  was 
next  increased  to  .2,  while  ^2  unchanged  at  2.  This  choice  parameters  produced  even  more 
severe  oscillations  in  the  predicted  backface  velocity  of  the  piatc.  At  this  point,  a  final  run  was 
made  using  the  PR0NTO2D  default  values  hj  =  0.06  and  hj  =  1.2.  Figure  3.3.4  describes  the 
results  for  model  A  (b^  =  0.06  and  1^2  =  1'2)  and  model  B  which  were  obtained  from  Figure 
3.3.1.  The  oscillations  in  the  time  history  for  the  fine  mesh  ,  tnodel  A,  have  completely  disap- 
peiued.  It  is  evident  that  the  later  time  results  for  the  fine  mesh,  model  A.  are  sensitive  to  the  aiti- 
cial  viscosity  parameters  and  hj-  Reference  [23]  contains  a  more  thorough  discussion  of 
artificial  viscosity.  However  at  tliis  time,  we  do  not  fully  understand  why  the  numerical  oscillations 
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occurred  for  higher  values  of  artificial  viscosity. 

The  future  inclusion  of  viscoelasticity/dispersion  (see  (2)  in  section  6,  Future  Directions)  in 
the  2D  and  3D  continuum  damage  tlieories  will  provide  additional  damping  that  should  help  pre¬ 
vent  such  oscillations  as  in  Figures  3.3.2  and  3.3.3.  We  note  that  the  WONDY  results,  with  up  to 
1500  zones  in  the  composite  plate,  included  dispeision  in  the  formulation  and  this  along  with  arti¬ 
ficial  viscosity  helped  to  prevent  the  oscillations  seen  in  some  of  the  PR0NT02D  models. 

Overall,  we  conclude  in  this  section  that  viscoelasticity/dispersion  may  have  an  important  in¬ 
fluence  upon  both  the  final  damage  profile,  and  in  the  capturing  of  the  shock  wave  as  it  propagates 
thru  the  thickness  of  the  composite  plates.  The  approach  discussed  in  [4]  will  be  extended  to  2D 
and  3D  situations.  With  regards  to  the  number  of  elements  thru  the  thickness,  it  appears  that  a  min¬ 
imum  of  100  may  be  required,  but  hundreds  of  elements  are  probably  not  necessary  to  achieve  suf¬ 
ficient  accuracy.  Certainly  more  analyses  still  need  to  be  performed  before  more  definitive 
conclusions  can  be  made.  However,  the  results  in  this  section  are  none  the  less  quite  encouraging. 


19 


PMMA  FLYER  -  intial  vel  s=  93nt/s 


«jt=0  all  nodes  Model  A  - 100  elements 

Model  B  -  25  elements 


7- 


CONTAa'  SURFACES  (x=0) 

_ ^  A- 


u^=0  all  nodes 


GRAPHn'E/PEEK  PLATE 
Model  A  -  400  elements 


Model  B  ■  100  elements 


FIGURE  3.1  -  PR0NT02D  model  of  ID  impact  of  PMMA  flyer  and 
graphite/peek  plate  (specimen  no.  204  in  16],  see 
Figure  2.1  also). 


TABLE  3.1  -  Material  and  damage  variables  used  for  the  PMMA  flyer  and  tlie 
grapliite/peek  plate^  in  tlie  PRONTO  analyses  -  SI  units. 


PMMA.(£lg;iigla2S).FlYfir 

E  (modulus)  =  9.595  XlO^Pa 

V  (Poisson’s  ratio)  =  .396 
p  (density)  =1185  kg/in^ 

Equation  of  State  -  Linear  Uj-  Up  Hugoniot  form  [8] 

Us*Co+s  Up 

Us=  shock  wave  velocity.  Up  =  particle  velocity 
Cq  =  2590.  s  =  1 .52  r  (Gmneisen  ratio)  =  .97 


GRAPHITE/PEEKPLATE 


Material  Properties 
p  (density) 

tnlanar  modulus^ 

sa 

£^,  (thru  tliickness  modulus) 
G®2  (transverse  sliear  modulus) 
V|2  (transverse  plane) 

(inplane) 

Damage  Parameters 


=  1579  kg/m^ 

=  69.X10^Pa 
=  13.45  X  10’  Pa 
=  7.X10’Pa 
=  .04 
=  .30 


a,  =  ctj  =  =  0.75,  =  0.95 

ilj  =  2.  X  10"^  «j  =  2 

112  =  11^=  1-4  X  10■^  ^2  =  "i  =  * 

a<.,o  =  70.  X  10^’  Pa  =  70.  x  W^Pa 

^C!0  ” 

^GsO  =  ^  ^ 

<Pc;-  )  “  ^G51  " 

Ultimate  strain  allowed 


^ULT 


=  0.015 
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DAMAGE 


FIGURE  3.2.1  Model  A,  fine  PR0NT02D  mesh  results. 


FIGURE  3.2.2  Mode!  B,  coarse  PR0NT02D  mesh  results. 


time  in  microseconds 


FIGURE  3.3.1  Predicted  backface  velocity  of  the  plate  for  model  1), 
coarse  mesh 


FIGURE  3.3.2  Predicted  backface  velocity  of  the  plate  for  model  A, 
Giic  mesh 
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— Model  A  -  400  elements  thru  thickness,  b2=2. 
— -  Model  B  <  100  elements  thru  thickness,  bj-.!,  b2=2. 


Time  in  microseconds 

FIGURE  33.3  Predicted  backface  velocity  of  the  plate  for  models  A  and 
B.  Note  oscillations  in  model  A. 


-  Model  A  •  400  elements  thru  tliickness,  bi3:.06,  b2=l>-2 

-  Modei  B  •  100  elements  thru  thickness,  bjs.l,  b2=2. 

100 

80 

Backface  Velocity 

of  Plate  ^0 

20 

I 

0* 

10  12  14  16  18  20 

Time  in  microseconds 

FIGURE  33.4  Predicted  backface  velocity  of  the  plate  for  models  A  and 
B.  No  oscillations  in  model  A. 
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3.3  Coarsf!  model  of  2D  impact  problem  -  circular  PMMA  flyer  and  graphite/ 
peek  plate 

Figure  3.4  describes  the  impact  experiments  that  were  conducted  in  [6],  which  consisted  of 
axisynimetric  PMMA  (plexiglass)  flyers  impacting  axisymmetric  composite  plates  composed  of 
either  graphite/peek  or  graphite/epoxy.  One  dimensional  models  of  experiment  204  are  discussed 
in  previous  sections.  See  Table  2.1  and  Figure  2.1  for  WONDY  model  results,  and  Figures  3.3. 1 
and  3.3.4  for  tlie  PRONT02D  one  dimensional  simulation. 

In  this  section,  we  will  discuss  the  2D  modeling  of  the  impact  experiment  using  PR0NT02D 
and  the  2D  continuum  damage  model  developed  in  [1-3].  In  general,  true  2D  effects  such  as  bend¬ 
ing  would  become  more  prevailent  at  later  time  during  the  experiment  and  simulation.  Early  time 
predictions  for  the  center  of  the  plate  would  be  dominated  by  one  dimensional  response  -  thus  the 
use  of  one  dimensional  WONDY  and  PR0NT02D  models  in  section  3.2.  Figure  3.5  is  a  listing  of 
the  formatted  PR0NT02D  input.  (Note  nodal  coordinates,  element  connectivity  and  boundary 
condition  data  are  given  in  another  input  file  and  is  unformatted  for  PR0NT02D).  Material  1  in 
Figure  3.S  represents  the  graphite/peek  plate,  while  material  2  is  the  PMMA  flyer.  The  material 
and  damage  variables  [1-3]  for  the  graphite/peek  plate  are  listed  in  Figure  3.5,  and  are  the  same  as 
contained  in  Table  3.1.  A  very  coarse  PR0NT02D  two  dimensional  model  is  indicated  in  Figure 
3.6  in  which  10  to  20  elements  have  been  used  thru  the  thickness  of  the  flyer  and  composite  plate. 
The  PMMA  flyer  contains  a  total  of 400  elements,  while  the  graphite/peek  composite  plate  is  dis¬ 
cretized  with  800  elements. 

In  Figure  3.7,  the  predicted  backface  velocity  of  the  centerline  of  the  plate  is  plotted  as  a  func¬ 
tion  of  time.  The  experimentally  measured  plate  backface  '/elocity  is  given  in  Figure  2. 1 ,  along 
with  WONDY  predictions.  In  addition,  PRONT02D  results  for  a  ID  simulation  using  100  and 
400  elements  thru  the  thickness  of  the  plate  are  given  in  Figures  3.3.1  and  3.3.4.  As  one  might  ex¬ 
pect  with  such  a  coarse  2D  discretization,  the  PRONT02D  results  of  Figure  3.7  represent  only  a 
crude  approximation.  Figure  3.8  describes  the  PRONT02D  predicted  thru  thickness  damage  in  the 
graphite/peek  plate.  Compare  this  result  to  Figure  2.1,  which  is  from  WONDY,  and  also  to  Figures 
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3.2.1  and  3.2.2  which  represent  PR0NT02D  one  dimensional  models  with  400  and  100  elements 
thru  the  thickness  of  the  plate.  In  general,  the  damage  results  in  Figure  3.8  are  veiy  approximate 
as  expected. 

Finally,  it  may  peihaps  appear  odd  to  refer  to  Figure  3.6,  in  which  1200  elements  (800  in  the 
plate)  have  been  employed,  as  a  coarse  mesh.  However,  one  is  trying  to  simulate  impact  experi¬ 
ments  in  very  hetergeneous  composite  materials.  Many  degrees  of  freedom  are  required  to  accu¬ 
rately  capture  the  shock  wave  as  it  progresses  thru  the  media.  Recall  in  section  3.2  model  B  (one 
dimensional  simulation)  with  1(X)  PRONT02D  elements  thru  the  thickness  compared  well  to  the 
early  time  experimental  results  for  the  backface  velocity  of  the  plate.  More  will  be  said  in  section 
5,  Summary  and  Future  Directions,  concerning  discretization  and  also  the  inclusion  of  viscoelas¬ 
ticity/dispersion  effects. 
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FIGURE  3.4  Impact  tests  conducted  in  [6]. 
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TITLE 

LAST  PROBLEM  -PMMA  ELVER, CR-PEEK  PL. -SHOT  3357-1200EL. 
PLANE  STRAIN 
BULK  VISCOSITY.. 1, 2. 

.ATERIAL.l, TRANS  ISO  W  DAM, 1579. 

PLAN.AR  MODULUS.  6  9.  E9 
THRU  MODULUS. 1 3 . 4SE9 
tIUTRANSPL..04 
NUINPLANE.  . 

SHEARTRANS.7.0E9 
Al.. 75 
A2..75 
A3. .75 
A4..9S 
ET1.2 . E-06 
111=2  . 

ET2.1.4E-0’ 

K2.1 

V1SIGG0<'.0.E$ 

VlTAUGOr.'O.EE 
VIPHIGO. . 5 
VLPHIG. . 5 
V2SIGG0.17O.Ee 
V2TAUG0.340.Ee 
V2PHIG0.). . 

VEPHIGd. 

EULT. . OIS 
END 

MATERIAL, 2, EP  HYDRODYNAMIC , 1185 . 

YOUNGS  HOOULUS.9 .S95E9 
POISONS  RATIO. .39fi 
YIELD  STRESS. 10. ElO 
H.TP.DElJISa  MODULUS.9 .  SSSEfl 
BETA.l . 

PRESSURE  CUTOFF. -1.E8 
END 

EQUATION  OF  STATE  2  .  MG  US-UP 
CO. 2590., S»l. 52, GAMMA.. 97 
END 

NO  DISPLACEMENT, X, 1 
INIT  VEL  MAT, 2,0.. -93. 

CONTACT  SURFACE,  100 , 200,  . 1 ,.  5 

TERMINATION  TIME.25.E-6 

OUTPUT  TIME, 1 . E-7 

PLOT  TIME, 1 .E-7 

PLOT  NODAL.VELOCITY 

PLOT  ELEMENT,  STRESS 

PLOT  STATE. VI . V2 , DOUT 

DEATH , 1 , DOUT , ABS ,  . 9  5 ,  . 1 

EXIT 

lEOBl 


FIGURE  3.5  PR0NT02D  input  for  coarse  mode!. 
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FLYER 


PLATE 


FLYER  -  .102  meters  radius,  .00345  meters  thick,  10X40  elements 

Plate  -  .102  meters  radius,  .02583  meters  thick,  20X40  elements 


FIGURE  3.6  Coarse  PR0NT02D  model  of  plexi-glass  flyer 
and  graphite/pcek  plate,  impact  vel.  93  in/s 
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Backface  Velocity 
of  Plate 


40 


20 


y  v  vvw 


5.  10 


0.00001  0.000015  0.00002  0.00002S 


TIME  IN  MICROSECONDS  (after  impact) 

FIGURE  3.7  Centerline,  backface  particle  velocity  for  the  coarse  PRONTO 
mesh  in  Fig.  3.6. 


Damage 


TIME  IN  MICROSECONDS  (after  impact) 


FIGURE  3.8  Thru  tliickness  centerline  damage  prediction  for  the  coarse 
PRONTO  mesh  in  Fig.  3.6. 


4.  DEVELOPMENT  OF  A  3D  CONTINUUM  DAMAGE  MODEL 
4.1  General  remarks 

The  development  of  a  robust  three  dimensional  (3D)  continuum  damage  model  for  transient 
analysis  of  thick  composites  has  proven  to  be  a  formidable  task.  A  simple  extension  to  tlie  2D  the¬ 
ory,  which  is  discussed  in  [1-3],  is  not  viable.  In  what  follows,  a  somewhat  different  approach  for 
3D  continuum  damage  is  developed. 

Essentially,  the  2D  continuum  damage  model  assumes  transverse  isotropy  in  the  1-2  plane  as 
indicated  in  Figure  4.1.  Delamination  damage,  representing  a  network  of  matrix  cracking  aligned 
in  the  1-2  plane  with  normals  in  the  3  direction,  is  denoted  by  VjCV,  in  [1-3])  and  depicted  in  Fig¬ 
ure  4.2.  Inplane  matrix  cracks,  which  may  traverse  or  lie  between  reinforcing  fibers  but  not  break 
fibers,  are  indicated  in  Figure  4.3  and  referred  to  as  damage  in  [1-3].  In  this  approach,  damage 
effects  were  tracked  and  applied  to  the  engineering  moduli,  shear  moduli,  and  Poisson’s  ratios. 
Specifically  for  the  coordinate  system  indicated  in  Figure  4.1,  the  following  damage  dependencies 
were  assumed: 

£33= 

~  ^22  ~ 

G31  = 

V31  =  (l-V2)(l-a3Vj)t,o, 

=  (l-a^vj)!)®^ 

where  the  superscript  “0”  denotes  undamaged  properties  and  the  fractions  0  <  a.  <  1 ,  i  =1,  4  are 
included  to  prevent  a  complete  loss  of  material  integrity  as  the  saturation  state  V  1  is  reached. 

The  assumption  of  transverse  isotropy  in  the  2D  fonnulation  represents  a  reasonable  approx¬ 
imation  that  allows  the  thni  thickness  direction  to  be  more  easily  discretized.  This  approximation 
eliminates  the  requirement  of  tracking  each  individual  ply  and  the  various  compliances  associated 
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with  tljem.  It  represents  essentially  a  homogenization  of  the  laminate.  Note,  however,  a  large  num¬ 
ber  of  elements  are  still  required  in  the  the  3  direction  to  accurately  model  transient  tliru  thickness 
effects  -  see  section  3.2. 

In  a  3D  model  of  a  laminated  composite,  the  transverse  isotropy  assumption  is  no  longer  valid. 
Total  damage,  which  will  be  represented  by  tlie  vector  D,  can  be  expressed  by: 

e  =  D,^r,■^D2e2  +  f>3^3  (4.2) 


where  g.  are  unit  base  vectors  and  are  individual  components  of  damage  with  normals  aligned 

in  the  f .  direction. 

A  direct  extension  of  the  2D  damage  model  to  3D  seems  quite  logical  at  first  glance. 
However,  this  approach  is  very  inefficient  because  it  would  necessitate  the  tracking  of  each  indi¬ 
vidual  ply  and  its  compliances. 


•  f2»  ♦v'  »a«4wcv  va«v  «4ia\a  i%  iv?  n»v 

stresses  [10-12]  as  opposed  to  the  compliances,  as  is  done  in  the  2D  model.  In  a  one  dimensional 
case,  the  effective  stress  a  with  damage  present  can  be  related  to  the  stress  of  the  undamaged  ma¬ 
terial  by; 


o  =  oxA/A  =  a/(l -D)  (4.3) 

where  D  =  ^  (4.4) 

A 

is  the  scalar  damage,  and  A  is  the  original  area,  and  A  is  the  effective  area  due  to  material  damage. 
From  eq.  (4.4),  D  is  interpreted  as  the  relative  reduction  in  area  caused  by  damage  due  to  micro¬ 
cracking.  For  general  anisotropic  damage  in  three  dimensions,  tnc  effective  stress  a  can  be  ex¬ 
pressed  in  a  tensorial  form  as: 


S  =  M{D)g  (4.5) 

where  M  is  the  second  order  material  damage  tensor  and  D  is  the  damage  vector  defined  by  eq, 
(4.2).  In  references  [10-12],  the  following  form  of  M  is  proposed  to  account  for  anisotropic  mate¬ 
rial  damage  in  the  principal  coordinate  system: 
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M  = 


(4.6) 

where,  as  in  the  2D  damage  theory,  the  fractions  0  <  <  1 ,  i  “1,2,3  have  been  added  to  prevent  a 

complete  loss  of  material  integrity  as  the  individual  components  of  damage  go  to  1. 

For  individual  plies,  eqs.  (4.5)  and  (4.6)  can  be  easily  applied,  since  the  principal  directions 
are  in  the  1  -2  plane  (see  Figure  4.1)  and  are  parallel  and  perpendicular  to  it.  Then  the  stress  tensor 
<j  can  be  transformed  to  the  global  1-2  axes.  In  this  work,  we  will  apply  eqs.  (4.5)  and  (4.6)  to 
groups  of  plies,  and  assume  the  additional  approximation  that  the  global  1-2  axes  represent  the 
principal  axes  (see  comments  at  the  end  of  this  section). 

For  linear  elasticity,  the  relationship  between  the  stress  and  the  strain  tensors  can  be  written  as: 

S  =  (4.7) 

wlicre  C  is  the  constitutive  tensor  and  e  is  the  strain.  In  the  case  of  an  orthotropic  lamina  in  which 
there  is  no  interaction  between  the  shear  stress,  eq.  (4.7)  is  of  the  form  [22]: 
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^12 

^13 
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^11 

^^22 

^21 

C22 

'^"23 

0 

0 

0 

^22 

^33 
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^32 

C33 

0 

0 

0 

^33 

«'l2 

0 

0 

0 

<^44 

0 

0 

®33 

<^13 

0 

0 

0 

0 

^^55 

0 

®13 

^23 

0 

0 

0 

0 

0 

^23 

®23 

(4.8) 


where  1,2,  and  3  are  the  principal  axes.  In  terms  of  the  engineering  constants,  the  components  of 
the  constitutive  tensor  Q  are: 

C|1  =  (1  -  V2jV22)/(E2^3p) 

^12  ”  (Vi2 “  ^21 

^13  ~  ('^IS  4'''i2''23^ “  ^31 

<^22  =  (l-v,3V3,)/(E,E3p)  (4.9) 

C23  -  (V23+V2,Vj3)/(£jE2p)  =  C32 

^33  “  (^  “''l2'^21^^ 

^44  ”  ^12'  ^55  “  ^13'  ^66  ~  ^23 

and 

P  =  (1-V,2V2,~V23V32-V3,V,3-2V2,V32V,3)/(£,£:2£3) 

where  v.^.,  and  represent  the  various  Poisson  ratios,  tlie  engineering  and  shear  moduli.  In 
addition,  there  is  also  the  following  relationship: 

v../£..  =  \./E.  (4.10) 

where  the  summation  convention  is  suspended. 
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Equations  (4.8)  and  (4.9)  can  be  constructed  for  eaclt  individual  ply.  However  as  discussed 

previously,  it  is  not  our  intent  to  track  tire  damage  in  each  ply  but  rather  in  groups  of  plies.  Tliis 

latter  approach  allows  more  flexibility  iit  the  discretization  of  tlie  tlim  thickness  direction  (3  axis). 

For  3D  and  2D  solid  finite  element,  explicit  transient  analysis,  such  as  tlie  programs  PRONT02D 

[8]  and  ABAQUS/EXPLICIT  [9],  employ  linear  solid  elements  with  reduced  or  one  point 

integration.  With  this  in  mind,  equations  (4.7)  -  (4.10)  can  be  applied  to  groups  of  plies  and  their 

effect  incorporated  into  a  single  element.  Within  a  linear  element  due  to  one  point  integration,  the 

stress  and  strain  do  not  vary  in  the  thru  thickness  direction  (or  in  the  other  two  directions).  Thus, 

eq.  (4.7)  for  a  group  of  N  plies  can  be  written  as: 

N 

(4.11) 

/■=  1 


In  lumping  together  groups  of  plies  with  orientations  a,b,c,d  for  instance,  we  appear  to  be  restrict¬ 
ing  the  number  of  element  allowed  in  the  3  direction.  As  discussed  in  sections  2  of  this  report,  in 
many  iiisiance.s  over  100  eleiucnis  thru  the  thickness  may  be  required  to  track  the  wave  nature  of 


the  response.  If  there  arc  only  10  groups  of  plies  with  orientations  of  a,b,c,d,  with  this  approach 
only  40  elements  in  the  3  direction  could  be  allowed  Ixjfore  individual  plies  would  have  to  be  mod¬ 
eled  with  more  tlian  one  element.  This  approach  would  be  very  difficult  particularly  with  respect 
to  pre-processing.  To  overcome  thus  defficiency,  an  additional  approx imatlion  is  introduced, 
which  assumes  each  group  of  plies  can  be  broken  down  into  subgroups  with  exactly  the  same  group 
of  plies  as  the  original  group  a,b,c,d,  but  with  scaled  thicknesses.  Thus,  if  ns  in  this  case,  there  are 
10  groups  of  plies  with  orientations  of  a,b,c,d  and  100  elements  thm  the  thickness  are  required,  it 
is  assumed  that  there  are  100  groups  of  plies  (a,b.c,d)  with  tluckness  1/10  of  the  original  layup. 
This  latter  homogenization  approximation  appears  to  be  a  reasonable  compromise  to  the  modeling 
of  each  individual  ply,  or  the  jx)ssibility  of  requiring  more  than  one  clement  per  ply  should  addi¬ 
tional  elements  be  necessary  in  the  3  direction. 


We  return  now  to  equations  (4.5)  and  (4.6).  The  strain  energy  W  of  an  undamaged  material 
can  be  expressed  as; 
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(4.12) 


Using  eq.  (4.7),  (4.12)  becomes: 

(4.130) 

in  tenns  of  the  strain  tensor  £  or 

w  =  iU3b) 

with  respect  to  tlie  stress.  The  strain  energy  Wp  for  a  damaged  material  is  expressed  as  [10]: 

Wp  =  (4.14) 

Using  cq.  (4.5)  produces: 

Wo  =  5£F^(M^rr'M)fif 
or 

h'd = y^c-'a 

where 

(4.15) 

Inverting  eq.  (4.15)  produces  the  damaged  constitutive  tensor: 

C  «  (4.16) 

Because  of  eq  (4.6),  M  is  assumed  to  be  a  diagonal  matrix,  and  thus  its  inverse  is  easy  to  obtain. 

Before  closing  this  section,  we  return  to  eqs.  (4.5)  and  (4.6).  The  M  matrix  is  assumed  to  be 
diagonal,  and  also  referenced  to  the  principal  coordinate  system.  In  general  for  laminated  compos¬ 
ites,  the  1  and  2  ditections  will  not  correspond  to  the  principal  directions  This  assumption,  how¬ 
ever,  can  be  partially  rectified  by  adjusting  the  thresholds  (see  section  4.2)  for  damage  to  more 
accurately  reflect  the  grouping  of  the  individual  plies  within  a  single  finite  element. 
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4,2  Details  of  the  damage  model 

Having  postulated  an  approach  to  model  3D  continuum  damage  in  a  thick  composite,  wc  next 
describe  the  specific  details  of  the  damage  model.  The  intended  application  is  for  composite  ma¬ 
terials  which  usually  contain  a  great  number  of  imperfection  sites  from  which  cracks  can  grow 
Thus,  the  nucleation  of  new  cracks  is  not  addressed.  Ratlier,  the  focus  is  on  the  growth  and  evolu¬ 
tion  of  damage. 

Recall  from  eq.  (4.2)  that  die  damage  Z?  is  assumed  to  be  a  vector  with  components  £>,  in  the 
1,  2,  and  3  directions.  As  in  the  2D  model  [1*3],  the  evolution  of  the  individual  components  of 
damage  are  assumed  to  be  governed  by  a  tlireshold  in  the  fonri: 


for  no  damage  growth 

(4.17) 

>0 

for  damage  growth 

where  (i.j)=1.2,3  and  there  is  no  sum  on  i:  F,  arc  scalar  threshold  functions;  <j  is  the  current  stress 
tensor;  and  is  an  array  of  cunent  tlireshold  parameters  which  are  functions  of  Dy,  the  damage 
component.  is  assumed,  as  in  [1-3],  to  be  of  the  Mohr-Coulomb  type  and  to  be  dependent  upon 

the  stress  measures  c,  and  t  (tension  and  shear)  as  follows; 

/  ^  2y/2 

=  i  +  (;^)  -^xr^yhi 

V  hi  J 

In  eq.  (4.18)  the  parameters  ai^e  related  to  specific  growth  threshold  strengths  and  t^.,  and  the 
Coulomb  friction  tangent  as: 

~  Ai~f2i 

(4.19) 

^Gi  ~  h/hi 

The  tension  and  shear  growth  thresholds  as  well  as  the  Coulomb  friction  tangent  in  eqs.  (4. 19)  are 
assumed  to  be  functions  of  the  damage  D-.  'I'he  threshold  surface  (F  -  0)defined  by  cq.  (4.18)  is 
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shown  in  Figui-e  4.4  along  with  the  various  parameters.  Note  that  in  this  figure  d  is  the  shortest 
distance  fron  an  external  state  (a,,  x)  to  the  threshold  surface. 


Using  eqs  (4.18)  and  (4.19),  we  further  develop  the  details  for  the  components  of  damage  D, 
and  its  evolution.  For  type  of  damage,  which  corresponds  to  delamination,  the  stress  measures. 
<5^.  and  x^.,  are  postulated  to  be  of  tlie  form: 

^f3  “  ®33 

X3  =  (<T^1  +  Ct22)‘^^  (4.20) 

The  grov/th  threshold  stresses  and  the  Coulomb  friction  tangent  in  eqs.  (4.19)  aie  then  taken  to  be: 

^G3  ”  ( ^  “  ^3)  ^G30 

1^3=  ('>■2') 

*Pg3  “  *Pc30'*'2^3(^C3I~^C30^ 

where  the  “0”  subscript  denotes  the  initial  undamaged  properties.  Note  that  for  0*^3  and  all 
resistance  to  damage  is  lost  as  D3  goes  to  1;  and  the  Coulomb  friction  tangent  tp^j  can  vary  lin¬ 
early  with  D]  as  die  damage  progresses. 


From  eqs.  (4.18),  (4.19),  the  variables  /23’ ^13  and  ^33  are: 


f23=V 


''G3 


- O 


- V^Gl^c3 

/l3  ~  ^C3‘^-^23 
/33  ~  f2^G3 


G3 


(4.22) 


To  complete  the  continuum  damage  formulation  for  JO3  damage,  evolution  equations  are  re¬ 
quired.  From  [1-2],  the  following  relationship  is  postulated  for  D3  (delamination)  damage: 

dD.  _ 

^  =  F-i{d^,D^)  (4.23) 

where  d^  is  the  shortest  distance  in  the  a,,  x  stress  space  (eq.  (4.20))  -  see  Figure  4.4  for  the  thresh- 
old  surface  F  =  0.  Forri3  -  0,  F  =  Oand^-  =  0  for  all  stress  points  interior  or  on  the  thresh- 
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old  surface  (F  =  0).  The  specific  form  of  eq.  (4.23)  employed  in  [1-2]  and  which  will  be  used 
here  is  expressed  as: 


(-is/Ocjo)  "■/<%(' -Oj)) 


(4.24) 


where  is  a  positive  power  exponent  for  (a  dimensionless  stress  distance),  iij  is  a  time 

constant  governing  the  rate,  and  the  term  1/(1-  d|)  causes  an  acceleration  of  the  damage  to 
complete  delamination  as  D3  goes  to  1 . 


Next  we  consider  Dj  and  D2  type  of  damage,  which  relate  to  matrix  cracking.  The  discussion 
will  focus  on  j ,  but  Dj  type  of  damage  will  be  assumed  to  be  of  the  same  form.  The  substituion 
of  2  for  the  index  1  in  the  equations  to  follow  will  yield  the  D2  formulation.  For  D^  type  of  dam¬ 
age,  the  stress  measures  and  Tj  required  foreqs.  (4.18)  and  (4.19)  are  postulated  to  be  of  the 
form; 


=  <^11 

Tj  =  +  (4.25) 

Equations  (4.25)  are  used  in  conjunction  with  eqs.  (4.18)  and  (4.19).  The  form  of  the  growth 
threshold  stresses  is  not  of  the  same  form  as  D3  damage,  W'hich  is  given  by  eq.  (4.19).  Dj  (and 
also  Dj)  damage  consists  of  an  ever  denser  network  of  cracks  which  tend  to  saturate.  There  is  no 
acceleration  to  a  catastrophic  failure  as  in  the  case  of  damage  which  represents  delamination. 
Saturation  can  be  forced  through  the  growth  thresholds  as  done  in  [1-2]; 

while  Coulomb  friction  tangent  is  taken  to  be 

‘Pgi  ~  ‘Pgio'*’ ^^Cii  “ (4.27) 

In  eq.  (4.26)  it  is  seen  that  as  D^  goes  to  1  damage  becomes  more  difficult  since  the  growth  thresh¬ 
olds  a^i  and  greatly  increase.  Reference  [1]  also  offers  another  form  for  the  threshold  func¬ 
tions  -  see  that  discussion  for  further  details. 
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For  D  j  (and  D2)  damage,  the  final  ingredient  required  is  the  specific  form  of  the  rate  evolution 
equation  which  is  described  by  eq.  (4.23).  Since  there  is  no  acceleration  to  a  catastropic  failure,  as 
in  the  case  for  D3  damage,  a  viable  form  is  taken  to  be: 

dDy 

J,  =  (rf/Oco)  '/’ll 

where  as  in  cq.  (4.24)  t)  ^  is  a  time  constant,  J,  is  normalized  with  the  virgin  growth  threshold 
stress  and  n  j  is  a  positive  term  for  the  dimensionless  stress  distance. 

Another  mode  of  continuum  damage  is  the  breakage  of  the  reinforcing  fibers.  To  account  for 
this  effect,  a  simple  maximum  strain  criteria  is  used  in  the  3D  model.  The  D^,  D2  and  Dj  forms 
of  damage  cause  softening  of  the  composite  response  and  thus  directly  influence  this  maximum 
strain  criteria. 

In  closing  this  section,  it  is  noted  that  at  the  time  this  repiort  is  being  written,  the  3D  continuum 
damage  model  is  being  programmed  as  a  user  written  subroutine  for  ABAQU S/EXPLICIT. 
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3 


figure  4.i. 


Portion  of  a  0  deg,  +  and  -  60  deg  ply  layup. 


figure  4.2  V3  dclamination  damage. 
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FIGURE  4*3  Vg  inplane  matrix  damage. 


FIGURE  4.4  Threshold  surface  for  the  onset  of  damage. 


5.  SUMMARY  AND  CONCLUSIONS 

The  results  of  this  report  document  progress  made  in  predicting  the  early  time  response  of 
thick  composite  plates  subjected  to  multi'dimensional  shock  loading.  In  general,  numerical  pre¬ 
dictions  compare  very  well  to  available  experimental  data.  Overall,  the  results  of  this  report  are 
very  encouraging  for  the  ID,  2D  and  3D  continuum  damage  theories  developed  at  NRL  for  the  ev¬ 
olution  of  damage  in  thick  composite  shell  structures. 

Dispersion  effects  in  a  ID  damage  model  are  discussed  in  section  2.  An  approach  is  developed 
that  employs  higher  order  derivatives  in  the  constitutive  relations.  This  is  shown  to  be  equivalent 
to  a  viscoelastic  theory.  This  approach  has  been  programmed  into  the  one  dimensional  finite  dif¬ 
ference  program  WONDY.  Comparisons  to  experimental  results  were  good.  However,  this  ap¬ 
proach  required  too  much  curve  fitting  of  data  from  individual  experiments  and  was  not  robust 
enough.  The  use  of  a  more  straightforward  viscoelasticty  theory  [15]  in  WONDY  along  with  the 
ID  damage  theoiy  proved  to  be  superior  and  is  documented  in  |4]. 

In  section  2.3,  using  the  WONDY  program  with  dispersion  and  damage  included,  ID  damage 
predictions  are  made  for  thick  composite  plates  subjected  to  underwater  shock.  Based  on  these  re¬ 
sults,  it  appears  that  spallation  is  not  a  real  threat  except  for  very  close  in  explosions. 

In  section  3,  PR0NT02D  finite  element  models  containing  the  2D  transversely  isotropic  con¬ 
tinuum  damage  theory  are  compared  to  experimental  and  WONDY  results  for  the  impact  of  plex¬ 
iglass  flyers  onto  graphite/peck  plates.  Plate  backface  velocity  predictions  from  the  ID 
PRONT02D  models  compared  veiy  well  to  experimental  measurements  [6]  and  very  finely  dis¬ 
cretized  WONDY  models  which  also  contained  viscoelasticty/dispcrsion.  One  hundred  4  node 
solid  elements  appear  to  be  sufficient  to  capture  the  magnitude  and  the  general  shape  of  the  shock 
wave.  The  inclusion  into  the  2D  continuum  damage  theory  of  viscoelasticty  as  used  in  [4]  as  well 
as  nonlinear  elasticity  in  the  thru  thickness  direction,  should  improve  these  results  further  -  see  sec¬ 
tion  6  item  (2)  for  additional  comments. 

A  new  3D  continuum  damage  theory  is  developed  in  section  4  for  thick  laminated  composite 
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plates.  The  3D  theory  represents,  in  part,  an  extension  of  the  2D  transversely  isotropic  damage  the¬ 
ory  contained  in  [1-3].  However,  three  dimensional  considerations  as  well  as  the  inefficiencies  of 
modeling  individual  plies  in  a  thick  composite  rcQuircd  a  slightly  different  approach.  A  3D  for¬ 
mulation  is  developed  which  applies  die  damage  vector  D  diiectly  to  the  stresses,  rather  than  the 
compliances  as  in  the  2D  theory.  In  general,  this  new  3D  damage  theory  appears  to  be  quite  prom¬ 
ising,  offering  real  potential  for  more  efficient  and  accurate  modeling  of  thick  laminated  composite 
plates  subjected  to  very  transient  loadings,  such  as  impact  and  underwater  shock. 
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6.  FUTURE  DIRECTIONS 

Based  on  the  success  of  our  efforts  at  NRL  to  date  concerning  the  development  of  numerical 
capabilities  for  predicting  the  response  of  thick  composite  shells  subjected  to  very  transient  load¬ 
ings,  we  plan  to  pursue  the  following  directions. 

(1)  The  2D  transversely  isotropic  continuum  damage  theory  [1-3]  will  be  implemented  into 
AB  AQUS/EXPLICrr.  This  capability  will  be  used  to  model  several  of  the  experimental 
impact  problems  discussed  in  [6]  in  which  plexiglass  flyers  were  fired  at  thick  graphite/ 
peek  plates  and  graphite/epoxy  plates.  Up  to  100  elements  thru  the  thickness  will  be  em¬ 
ployed.  True  2D  problems  in  which  the  diameter  of  the  flyer  is  smaller  tl^'m  the  plate 
(edge  effects)  as  well  as  bending  and  transverse  shear  effects  (later  time)  will  be  investi¬ 
gated.  Direct  comparisons  to  the  experimental  data  in  [6]  will  be  made,  especially  plate 
backface  velocity  histories.  In  addition,  damage  predictions  will  be  compared  to  any 
available  data  such  as  contained  in  [5]. 

(2)  Viscoelasticity  is  important  in  the  response  of  composite  plates  to  dynamic  loading. 
Therefore,  a  viscoelastic  formulation  similar  to  [15,4]  will  be  incorporated  into  the  2D 
continuum  damage  theory  and  implemented  into  ABAQUS/EXPLICIT.  This  formula¬ 
tion  will  of  course  also  introduce  dispersion  into  the  composite.  Also  included  will  be 
nonlinear  elasticity  in  the  normal  or  thru  thickness  direction  of  the  plate. 

(3)  The  new  3D  continuum  damage  theory,  developed  in  section  4  of  this  report,  will  be  fully 
implemented  into  ABAQUS/EXPLICIT.  Viscoelasticity  and  nonlinear  elasticity  in  the 
thru  thickness  direction  of  the  plate  will  also  be  introduced, 

(4)  A  nev/  thick  plate/shell  continuum  damage  theory  based  upon  the  new  3D  theory  of  sec¬ 
tion  4  will  be  developed. 

(5)  We  will  begin  to  investigate  the  introduction  of  the  2D  and  3D  continuum  damage  theo¬ 
ries  into  an  implicit  FEM  program  such  as  ABAQUS/STANDARD.  Implicit  codes  have 
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typically  been  much  more  useful  than  explicit  programs  for  structural  dynamics  types  of 
problems,  such  as  underwater  shock,  in  wliich  one  is  not  trying  to  track  the  actual  wave 
propagation  in  the  structure.  Projectile  impact  types  of  problems,  on  the  other  hand,  aie 
an  example  best  modeled  usually  with  an  explicit  program.  Implicit  approaches  require 
the  fonnulation,  computation  and  inverse  of  a  global  stiffness  matrix,  but  can  use  much 
larger  time  steps  and  are  usually  unconditionally  stable.  The  integration  of  the  governing 
constitutive  relations  is  also  usually  more  involved  in  implicit  formulations. 
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8.  APPENDICES 


8.1  Appendix  A  >  Reduced  ID  damage  model  implemented  in  WONDY 

As  mentioned  in  the  Introduction.  WONDY  has  been  used  extensively  in  our  work  to  study 
ID  plane  strain  behavior.  In  this  appendix,  the  details  of  the  ID  continuum  damage  model  which 
has  been  programmed  into  WONDY  are  given.  This  constitutive  model  stems  from  the  2D  damage 
theory  developed  in  [1-3],  and  also  includes  the  nonlinear  response  of  the  material  in  compression. 
Typical  material  and  damage  parameters  are  given  for  graphite/peek  composite  plates  [6]. 

In  the  one  dimensional  WONDY  computer  models,  the  1  direction  is  thru-thc*thickness  of  the 
composite  plate,  while  the  2  direction  is  inplane  (ID  plane  strain  models).  Tire  constitutive  equa¬ 
tion  is  given  as  follows: 

Ojj  =  Cn(V,)ejj  in  tension  (8. In) 

Ojj  =  JSfiEjj  + /(r.2ej, -t- 1(^36^ j  in  compression  (o.ll)) 

where  <jjj  is  the  stress,  £,j  is  the  strain,  Cj,  is  the  material  stiffness,  V^  is  tlte  dmnage  thru  the 
thickness  (delantination)  and  the  K*s  are  constants  interpolated  from  the  data  in  [6].  Typical  val¬ 
ues  used  for  the  K's  for  graphite/peek  composite  plates  ore; 

=  Cl,  (0)  =  11.7  GPa 

K2  -  114.  GPa  (8.2) 

A3  =  684.  GPa 

The  damage  model  based  upon  [1-3]  is  expressed  as: 

(1-V23)E„(V,) 

“  ‘  ( 1  -  V23  -  2x^2  ( V,)  (E22/  (C,  1  ( V, ) ) ) ) 

£,,(V,)  =  (1-V2)£«,  (8.3) 

v,2(V,)  =  il-V])vl2 
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in  wliich  Ej,  and  £22  engineering  moduli,  and  v,2  and  are  the  Poisson’s  ratios.  The 

evolution  of  the  damage  V,  is  in  the  fonn  [l-SJ: 


£V, 

(it 


^GO 


"i 

/(Tljd-Vj)) 


di(V,,/)  -  WaA[0,Cj,-o^] 


(8.4) 


In  eqs  (8.4),  t  is  the  time  in  seconds  and  d,,  o^q,  cr^  arc  defined  in  [1-3],  Typical 

values  used  in  eqs.  (8.3)  and  (8.4)  for  graphite/peek  are; 


£?i  =  13.45  GFa 
E22  ~  69.  GPa 
v^2  =  0.04 
V23  =  0.30 

and  for  the  damage  threshold  parameters 

^GO  “ 

rtj  =  2. 

1),  =  4.  X  10"^  see 


(8.5) 


(8.6) 
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8.2  Appendix  B  -  Listing  of  the  ABAQUS/EXPLICIT  user  written  subroutine 
for  2D  continuum  damage. 
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SUBROUTINE  VUMAT ( 

C...Read  only  variables 

1  NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL, 

2  STEP_TIME,TOTAL_TIME,DT,CMNAME, PROPS, DENSITY , 

3  STRAIN_INC , REL_SPIN_INC , 

4  TEMP_.OLD ,  STRETCH_OLD ,  FI  ELD_OLD , 

5  STRESS_OLD , STATE_OLD , ENER_INTERN_OLD , ENER_INELAS_OLD . 

6  TEMP_NEW,STRETCH_NEW,FIELD_NEW, 

C... Write  only  modifiable  variables 

7  STRESS_NEW ,  STATE_NEW ,  ENER_INTERN_NEW ,  ENER_INELAS_NEW ) 


C 


INCI/JDE  '  ABA_PARAMS .  INC ' 


C...Thi3  rountine  is  the  2d  continuum  damage  theory  developed 
C  by  Nemes  and  Randles  and  in  the  PRONt62D  code  - 
C  Transversely  Isotropic  Damage  Model.  Note  as  in  PRONT02D 

C  the  1,2,3  axes  in  ABAQUS  are  not  the  same  as  given  in  the 

C  report  and  paper  l->2,  2->3,  3->l  (from  paper  to  ABAQUS 

C  and  PRONTO) .  Nomenclature  very  similar  to  what  Nemes  used 

C  in  PR0NT02D  has  been  employed.  C.T.  DYKA  10/92 
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 


...The  7  state  variables  are  stored  as  for  element  I: 
applies  for  old  and  state  variables 

STATE_OLD ( 1 , 1 )  =vl  (Thru  thic)cness  damage) 

STATE_OLD ( I , 2 ) =V2  (In  plane  damage) 

STATE_OLD ( I , 3 ) =DOUT  (Switch  for  the  death  option) 

.  .  . 

...NOTE;  Not  sure  ABAQUS /EXPLICIT  ha...  regular  DEATH  OPTION!!!! 

ABAQUS  does  not  track  the  strains,  only  the  increments,  so 
use  the  state  variables  to  track  the  total  strains  -  which 
we  need  to  determine  when  the  ultimate  strain  has  been 
reached. 


STATE_OLD(I,4)=total  strain (1,1) 

STATE_OLD(I, 5)=total  strain (I, 2) 

STATE_or,n  (1,6)  stotal  strain  (1,3) 

STATE_01jD  (1,7)  =total  Strain  (1,4) 
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 
C...A11  arrays  dimensioned  by  (*)  are  not  used  in  this  subroutine. 


DIMENSION  PROPS (NPROPS) , DENSITY (NBLOCK) , 

1  STRAIN_INC (NBLOCK, NDIR+NSHR) , REL_SPIN_INC (* ) , 

2  TEMP_OLD ( * ) , STRETCH_0LD ( * ) , 

3  FIELD_0LD(*) , 

4  STRESS_0LD (NBLOCK, NDIR+NSHR) , STATE_OLD (NBLOCK, NSTATEV)  , 

5  ENER_INTERN_OLD ( * ) , ENER_INELAS_OLD ( * ) , 

6  TEMP.  NEW ( * ) , STRETCH_NEW ( * ) , 

7  FIELD_NEVJ(*)  , 

8  STRESS._NEW  (NBLOCK ,  NDIR+NSHR)  ,  STATE_NEW  (NBLCCK,  NSTATEV)  , 

9  ENER_INTERN_NEW  (  *  )  ,  E!JER_IMELAS_NEW  (  * ) 

CHARACTER *8  CMNAME 

C###  ###########  ######  #()###«####((############  ######0  0 

C...  Material  Properties 


Elio  ::: 

PROPS ( 1 ) 

E22''  = 

PROPS ( 2 ) 

rnu.j:  0 

=  PROPS (3) 

RNU130 

=  PROPS (4) 

G210 

::  PROPS  (5) 

A1 

=  PROPS (6) 

A2 

=  PROPS (7) 
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A3  =  PROPS (8) 

A4  =  PROPS (9) 

ETAl  =  PROPS (10) 

ENl  =  PROPS (11) 

ETA2  -  =  PROPS (12) 

EN2  =  PROPS  (13) 

SIGGO  =  PROPS (14) 

TAUGO  =  PROPS (15) 

PHIGO  =  PROPS (16) 

PHIGl  =  PROPS (17) 

SIGGOB  =  PROPS (18) 

TAUGOB  =  PROPS (19) 

PHIGOB  =  PROPS (20) 

PHIGIB  =  PROPS (21) 

EULT  =  PROPS (22) 

DO  100  I=l,NBLOCK 

. . .  Continuum  damage  variable  and  DOUT  to  track  death 
. . .  of  an  element 
■  VI  =  STATE_OLD(I.l) 

V2  =  STATE_OLD(I,2) 

DOUT  =  STATE_OLD (1,3) 

.  . .OPTION  TO  KILL  AN  ELEMENT  WITHOUT  THE  DEATH  OPTION 
IF (DOUT  .GE.  .98)V1=.98 
IF (DOUT  .GE.  .98)V2=.98 


. . .  Stresses  (abreviations) 

SIGl  =  STRESS_OLD(I,l) 

SIG2  =  STRESS_OLD(I,2) 

SIG3  =  STRESS_OLD ( I , 3 ) 

SIG4  =  STRESS_OLD(I,4) 

-  UPDATE  ELASTIC  PARAMETERS.  Note  the  1,2,3  directions  are  switched 

from  the  reference  report  and  paper 
E22=E220* (l.-Vl**2) 

E11=E110* (1 . -A1*V2**2) 

G21=G210* (l.-Vl**2)*(l.-A2*V2**2) 

RNU21=RNU210* (1 . -Vl**2) * (1 . -A3*V2**2) 

RNU13=RNU130* (1-A4*V2**2) 

C22= (1.-RNU13)*E22/ (1 . -RNU13-2 . *RNU21**2*Ell/E22 ) 

C12 =RNU2 1 *E1 1 / ( 1 . -RNU13 -2 . *RNU2 1 * * 2 *E1 1 /E2 2 ) 

Cll= (1./ (1.+RNU13) ) *(1.-RNU21**2*E11/E22)*E11/ (1 . -RNU13-2 . * 

1  RNU21*'‘'2*E11/E22) 

C13= (1 . / (1 .+RNU13) ) * (RNU13+RNU21**2*E11/E22) *E11/ (1 .-RNU13-2 . * 

1  RNU21**2*E11/E22) 

C44=2.*G21 


- UPDATE  STRAINS  -  Note  we  use  the  state  variables 

STATE_NEW(I,4)  =  STATE_OLD ( I , 4 )  +  STRAIN_INC ( I , 1 ) 

STATE_NEW(I,5)  =  STATE_OLD ( I , 5 )  +  STRAIN_INC ( I , 2 ) 

STATE_NEW(I,6)  =  STATE_OLD ( I , 6 )  +  STRAIN_INC ( I , 3 ) 

STATE_NEW ( I , 7 )  =  STATE_OLD ( I , 7 )  +  STRAIN_INC ( I , 4 ) 

.  , .Determine  the  largest  strain  in  each  element  -for  the  DEATH  option. 
. ..  Look  at  strain (1,1)  and  strain (3, 3)  components  only 
STRC  =  MAX ( STATE_NEW (1,4), STATE_NEW (1,6)) 

-  COMPUTE  VI 

IF(V1  .GE.  .98)GO  TO  20 
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oo  no  no  o  nonon  n  non 


SIGG  -  SIG^! 

TAU  =  SIG4 

SIGGG  =  SIGGO*  (l.--Vl**2) 

TAUG  =  TAUGO* (l.-Vl**2) 

PHIG  =  PHIGO  -  (PHIGO-PHIGl) ♦Vl**2 

CALL  DCALC  (SIGGG. TAUG, PHIG, SIGG, TAU, Vl , D1 , N) 

IF(N.GT.O)  GO  TO  105 

IF  (Dl.LT. 0.00001)  GO  TO  10 

VID  =  (  (D1/SIGG0)**EN1)  /  (ETA1-*(1.-V1**2)) 

VI  =  VI  +  V1D*DT 
IF(Vl.GT.n.98)Vl=. 98 

-  COMPUTE  V2  (also  called  VS) 

]0  IF  (V2.GE.0.98)  GO  TO  20 
SIGG  =  .5* (SIG1+SIG3) 

TAU  =  SQRT(SIG4**2  +  .25*  (SIG1-SIG3*''2) 

SIGGB  =  SIGGOB/ (1. -V2**2) 

TAUB  =  TAUGOB/ (1-V2**2) 

PHIB  =  PHIGOB  -  (PHIG0B-PHIG1B)*V2**2 
CALL  DCALC (SIGGB, TAUB, PHIB, SIGG, TAU,  V2 , D2  ,  N) 

IF(N.GT.O)  GO  TO  105 
IF(D2 .LT. 0.00001)  GO  TC  20 
V2D=  (D2/3IGG0B)**EN2/ETA2 
V2  =  V2  +V2D'-DT 
IF  (V2.GE.0.98)V2=.98 

.  ..NEMES  ALSO  USED  -  I'LL  INCLUDE  V2  DAMAGE  TO  KILL  ELEMENT 

20  IF(V1  .GE.  .98  .OR.  STRC  .GT.  EULT  .OR  V2  .GE.  .98)DOUT=1.0 

20  IF ( STRC. GT. EULT) DOUT=l .0 
-  UPDATE  STRESSES 

STRESS._NEW  ( 1 , 1 )  =C1 1  * STATE_NEW  ( 1 , 4 )  +C1 2  * STATE_NEW  ( 1 , 5 )  + 

1  C13  *STATE_NEW  ( ,  6 ) 

STRESS_NEW  (1,2)  =C2  2 * STATE_NEW  ( 1 , 5 )  +C1 2  *  ( STATE_NEW  ( 1 , 4  )  + 

1  STATE_NEW;i, 6) ) 

STRESS_NEW ( 1 , 3 ) =C1 2  *  STATE_NEW  (1,5) +C1 3  *STATE_NEW  ( 1 , 4 )  + 

1  Cll *STATE_NEW (1,6) 

STRESS_NEW (1,4) ^C44*STATE_NEW (1,7) 

...Update  state  variables  -  damaae  Vl,V2  and  DOUT  (element  death) 
ST.\TE_NEW  ( 1 , 1 )  =  VI 
STATE_NEW (1,2)  =  V2 
STATE_NEW ( I , 3 )  =  DOUT 
100  CONTINUE 

RETURN 
105  CONTINUE 

...STOP  CALCULATIONS  N  IS  OUT  OF  BOUNDS 
STOP 
END 

SUBROUTINE  DCALC  ( SIGG, TAUG, PHIG, SIG, TAU, V,  D, N) 

.  . .  THIS  SUBROUTINE  DOES  THE  THRESHOLD  CALCULATIONS  FOR  VI  AND  VS 
(V2)  DAMAGE.  EACH  CALL  IS  FOR  ONE  OR  THE  OTHER 
D=0  . 

N=0 

SIGG=SIGG0* (1 . -V**2 ) 

TAUG=TAUG0* (1 . - V**2 ) 
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C  PHIG=PHIGO- (PHIGO-PHIGl) *V**2 
FO= (TAUG/PHIG) **2/SIGG 
F1=.5*(F0+SIGG) 

F2=:.5*(F0-SIGG) 

F3=PHIG*F2 

IF( (F1-SIG)/F2.GE.SQRT(1.+(TAU/F3)**2) >RETURN 
C1=F2* (F1~SIG)/(F2^*2+F3**2) 

C2=F3*TAU/ (F2**2+F3**2) 

IF(Cl.LT.O.)  GO  TO  10 
X=C2 

5  X1=C2+C1*X/SQRT(1.+X**2) 

N=N+1 

IF(N.GT.IOO)  RETURN 
XDIFF=ABS (XI -X) 

X=X1 

IF(XDIFF.GT.0.0001)GO  TO  5 
GO  TO  20 
10  X=C2/(1.-C1) 

15  X2= (l.+X**2) **1.5 

Xl=  (C2+C1*X’‘*3/X2)  /  (1.-C1/X2) 

N=N+1 

IF  (N.GT. 100) RETURN 
XDIFF=ABS(X1-X) 

X-Xl 

IF  (XDIFF.GT. 0.0001)  GO  TO  15 
20  TAOO  =  F3*X 

SIGO  =F1-F2*SQRT(1 .+X**2) 

D=SQRT( (SIG-SIGO) **2+ (TAU-TAUO) **2 ) 

N=0 

RETURN 

END 
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